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The Brusselator is a generic reaction-diffusion model for a tri-molecular chemical reaction. 
We consider the case when the input and output reactions are slow. In this limit, we show 
the existence of AT-periodic, spatially bi-stable structures, mesas, and study their stability. 
Using singular perturbation techniques, we find a threshold for the stability of K mesas. 
This threshold occurs in the regime where the exponentially small tails of the localized 
structures start to interact. By comparing our results with Turing analysis, we show that 
in the generic case, a Turing instability is followed by a slow coarsening process whereby 
logarithmically many mesas are annihilated before the system reaches a steady equilibrium 
state. We also study a "breather" -type instability of a mesa, which occurs due to a Hopf 
bifurcation. Full numerical simulations are shown to confirm the analytical results. 



1 Introduction 

In 1952, Turing proposed that the formation of spatial patterns during morphogenesis 
could be explained in terms of the instability of a homogeneous steady-state solution of 
a reaction-diffusion network describing the evolution of a set of morphogens [31]. Turing 
himself illustrated his ideas on two chemical models. Turing's original work is primarily 
concerned with the stability analysis of the homogeneous steady-state solution of the 
rate equations for the interacting morphogens [21]. The main point of biological interest, 
however, is whether stable spatial structures may be generated beyond the instability, 
i.e., whether the rate equations admit stable (and positive) inhomogeneous solutions 
exhibiting the most characteristic features of morphogenetic patterns. 

This point has been taken up seriously in the early seventies. More systematic numer- 
ical studies of Turing's model have been performed showing irregular spatial structures 
[4]. This led to serious reservations about the relevance of Turing's theory in develop- 
mental biology, particularly its ability to generate regular patterns. But these criticisms 
originate from a somewhat unfortunate choice of Turing's example and they do not touch 
the essential points of Turing's theory. More specifically, the obvious requirement that 
the rate equations must admit positive and bounded solutions is not satisfied in Turing's 
example [11]. 
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Figure 1. (a) Turing instability and mesa-type localized structures. The parameters are A = 
1, B = 8, e = 10 -4 , D = 10, t — 10. Initial condition was set to u = A, v = A/B, perturbed 
by a very small random noise. Note the logarithmic scale for time. Initial Turing instability 
triggers a k — 7 mode at time t ~ 50. Thereafter a coarsening process takes place until there 
are only two mesas left, (b) Slow-time oscillatory instability of a single spike solution to (1.4). 
The parameter values are A = 1, B = 8, D = 10, £ = 0.00025, r = 0.999. 



Because of this misunderstanding on Turing's fundamental contribution, other two- 
variable models satisfying the law of mass action have been explored. In 1968, Prigogine 
and Lefever [28] introduced a two variable system that exhibits an autocatalytic reaction 
(called the "Brusselator"). The simplicity of the rate equations motivated analytical and 
numerical studies which showed the existence of stable structures [3] . The Brusselator is 
based on the following intermediate reactions for the two chemical intermediates X and 
Y: 

A^X, C + X^Y + F, 2X + Y^3X, X -> E. (1.1) 
The global reaction is A + C — > F + E and corresponds to the transformation of inputs 
products A and C into output products F and E. We assume (without loss of generality) 
that the rate constants for the first and last step are equal to r whereas the intermediate 
rate constants are one. The rates equations then become 

-rA-CX + X 2 Y - rX 
CX ~ X 2 Y. 



X f — D T X, r 



(1.2) 
(1.3) 



In this paper we concentrate on the case where the input and the output reaction 
steps are slow in comparison to the intermediate steps, so that r is small. In the absence 
of diffusion terms, this implies that the total concentration X + Y is a slowly changing 
quantity. Linear stability reveals the presence of a periodic orbit which arise through a 
Hopf bifurcation of the steady state X = A,Y = C/A. In the limit r — > such a periodic 
orbit consists of a slow and a fast phase, so that a pulse train can be observed. 
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There is a large body of literature on the space-independent Brusselator (D x = = D y ) 
and its extensions. Under different assumptions on the parameters, the equilibrium state 
is only marginally stable and an addition of small random perturbations (or small periodic 
forcing) can lead regular to large-amplitude pulses [27], [20]. Other authors have looked 
at coupled brusselator systems which can exhibit chaos and syncrhonization [32], [40], 
[6]. 

Since the discovery of spatial patterns in 1970's, various Turing patterns in the Brus- 
selator were studied both numerically and analytically in one, two and three space di- 
mensions. These include spots, stripes, labyrinths and and hexagonal patterns [12], [21], 
[36], [30], [37], oscillatory instabilities and spatio-temporal chaos [38], [39]. While Turing 
analysis and its weakly nonlinear extensions have been successful at detecting and classi- 
fying possible pattern types, its range of applicability is limited. Indeed Turing patterns 
are assumed to be small sinusoidal perturbations of a homogeneous state. In practice 
however, many patterns are localized and contain sharp transitions such as spikes and 
kinks. As we will see, this is particularly the case under the assumption of small r. In 
this paper we study mesa-type patterns. 1 These are box-like patterns that join two flat 
regions of space with sharp transition layers, such as shown in Figure 2. Such patterns 
are not amenable to Turing analysis since they are far away from the homogeneous state. 

Previous mathematical efforts on the formation of localized patterns concentrated on 
two variable reaction-diffusion models where the dynamics is controlled by the interaction 
between a slow and a fast variable. For the Brusselator without diffusion terms and with 
r small, the total concentration X + Y plays the role of a slow variable. As a result, 
localized patterns patterns are expected when diffusion is permitted. 

For our analysis it will be convenient to use the following scaling, 

ru t — Deu xx + eA — Bu + u 2 v — eu 
v t = Dev xx + Bu - u 2 v, 



(1.4) 



where 



^ * = r^, B = C%L, D=^, (1.5) 

J_S X J_S X J_J X I 



u = X, v = tY (1.6) 
and the spatial domain is x £ [0,1] with the zero flux boundary conditions 

u x = = v x at x — and 1. (1.7) 

In this work we assume the following conditions, 

r < D x < 0(D y ) < 1, A = 0(1), C = o(^j. (1-8) 

In terms of our scaling we have 

e£><l, £>>1; 0(A) = 1 = 0(B), (1.9) 

0(t) > I. (1.10) 

1 Mesa means table in Spanish; it is also a name given to square boulders which are found 
in the Colorado desert. The use of this term was suggested by Fife [13]. 
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As a motivation, let us present two numerical examples. On Figure l.a we show a typical 
time evolution over a very long time interval. The homogeneous steady state is unstable 
because of Turing instability but the expected Turing sinusoidal pattern only appears 
after a delay (t > 50). Turing's structure then gradually deteriorates and relatively 
quickly moves into a new pattern formed by several localized mesa-type structures. They 
then undergo a coarsening process over a logarthmically long time-scale and, eventually, 
only two mesas remain. Turing's analysis can be used to predict the first pattern (t ~ 50) 
but it cannot anticipate the coarsening process or the number of final mesas. As we shall 
demonstrate, in the case B > A 2 , the Turing pattern is characterized by a wave number 
proportional to k = O while the number of stable long-time mesas has the order 

K = O ( ^oli n i ) ' so tnat ^ 1S logarithmically bigger than K. This is the underlying 
cause of the coarsening process observed in Figure l.a. 

In the second experiment shown in Figure l.b, a single mesa undergoes a "breather" - 
type oscillatory instability which eventually leads to its extinction. Both coarsening and 
the breather instability occur at a slow timescale. The main goal of this work is to describe 
these instabilities analytically. 

Our results are related to the study of bistable systems, see for example [29], [14], [22] 
[23], [24], [25], [26], [8]. Mesa patterns also appear in the FitzHugh-Nagumo model [14], 
certain phase separation models such as Cahn-Hilliard, Allen-Calm, [1], [2], [7] and block- 
copolymers [29]. For these systems the resulting spectral problem has small eigenvalues, 
also called critical spectra, that tend to zero with the thickness of the interfaces. Typically 
k such layers are stable [25]. However as we show in this paper, if the number of layers 
is excessively large, instabilities can occur. This happens when the exponentially small 
interaction outside the interface locations cannot be ignored. Our main new contribution 
is to study this interaction, and to show that it has a destabilizing effect. 

1.1 Summary of main results 

We now summarize our main results. We first describe the shape of the equilibrium K 
mesa solutions. In Section 2 we show the following. 

Proposition 1 Consider the equilibrium- state problem, 

= sDv xx + Bu - u 2 v, = sDu xx + sA + u 2 v - (B + e) u, x e [0, 1] , (1.11) 
u' = 0=v' at x = 0andx = l. (1.12) 

in the limit (1.9) and suppose that 

A 2 < 2B. (1.13) 

Then there exists a K-mesa symmetric solution to (v,u) of the following form. 
Let 

w = 3^/B/2, l = J^^, d= J(- l > (L14) 

x u =Xi- 1 -, x ri = Xi + ^, x l = - 2 ^ ^ for i = l...K. (1-15) 




0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 



Figure 2. An example of a three-mesa equilibrium solution for v. Here, 
K = 3, A = 2, B — 18, eL> = 0.02 2 , w = 9, Z = 0.11. 



For x away /rom x r i , , we /lave: 



V ~ < 

I 3 



.t e [0, 1]\ U [x/^av 

— . x £ U (xn 7 X T i^ 



For x near the interfaces x ril xu we have, 

vi(x-xu), x-xu<0(^§) 

Vj-i^pC •fori ) i *^ — *^v% ^ ^(\/ ) 



w/iere 

2 1 
u r = W7 - + w - tanh 

Finally, 



Wq x 



wq x 



2 1 , / 

I . vi = wo - — wo~ tanh 

3 V^dJ 3 °3 V 3 



It ~ — f • 



(1.16) 
(1.17) 

(1.18) 
(1.19) 



A typical such solution with K = 3 is illustrated on Figure 2. 

We next analyse the stability of such equilibrium. There are two distinguished limits of 
interest, either DK 2 <C O ( ln a e ) or DK 2 = O ( ^ - ) . The former is studied in Sections 
3 and 4 while the latter in Section 5. In Section 3 we derive rather precise results for 
eigenvalues, summarized in the following theorem. 

Theorem 2 Consider a K mesa solution of Proposition 1. Suppose in addition that 

1«M 2 «0(4-] and O (r - 1) > 0. (1.20) 

V em e J 
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Such solution is stable when t-1 > and unstable when t — 1 <C 0. There are 2K small 
eigenvalues of order O (e) ; all other eigenvalues are negative and have order > O (De) . 
The smallest 2K eigenvalues are given by 



-1 ± J I - 2K 2 dl fl - cos ( 2£)1 

A ? ± ¥ ; — e /or j = l.../T-l; (1.21) 

3± 2 (r - 1) 3 J . v ; 

—if? -1 

A_ - e , A+ = -e. (1.22) 

t — 1 r — 1 

one? are aiZ negative when r > 1 , and positive when t < 1 . TTie transition from stability to 
instability occurs via a Hopf bifurcation as t is decreased past Th where to leading order, 

7fc~l. 

Note that near the interfaces, the gradient changes on the order 8 — \J eD. In terms of 
8, the scaling DK 2 — 0( el *a £ ) can be written as 

D = 0(^exp(^)). (1.23) 

Thus Theorem 2 confirms the stability of K mesas when r > 1 as long as D is not 
exponentially large in 8. In the contrary case, we derive the following result in Section 5. 



Theorem 3 Suppose that 



and let 



t > 1 



(1.24) 



where D\ <~ < 



2 / 12v / 2AB 3 / 2 



2e In' 



(V2B-A)' 

(V2B-A) 2 
2s ln 2 ( 1^B3/A 



2A 2 < B 



2A 2 > B 



+ l.s.t. 



(1.25) 



Here, l.s.t. denotes logarithmically small terms. Then a K mesa symmetric equilibrium 
with K > 2 is stable if D < Dk and is unstable otherwise. Moreover, a single-mesa 
equilibria K = 1 is always stable. A more precise value for D\ is given in Proposition 8. 



Theorem 3 states that the instability threshold occurs when D is exponentially large. 
In this case the exponentially fast decay outside the interface locations must be taken 
into account. Their exponentially weak interaction is responsible for an eventual loss of 
stability. Figure 3 illustrates this proposition. Indeed using the parameters used in that 
simulation we deduce from Proposition 8 that D\ = 20.96; so that D 2 = 5.28. Since 
D = 10 > D2, the two-mesa equilibrium state is unstable. 

Our next result is about the presence of a Hopf bifurcation when r is near 1. 



Theorem 4 Suppose that 



< DK 2 < O 



e In 2 £ 



(1.26) 
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Figure 3. Slow-time competition instability of a two-mesa solution to (1.4). The parameter 
values are A = l, 5 = 8, D = 10, e = 0.00025, r = 1.3 



Let 



(1.27) 



where K,d,l are as in Proposition 1. Then a K-mesa solution undergoes a Hopf bifurca- 
tion when t — T) l+ . It is stable when r > t%, and unstable otherwise. When r = Th + , the 
corresponding eigenvalue has value 



A+ - iV8K (e 3 DB) 



1/4 



(1.28) 



Figure l.b illustrates the type of oscillations that occur when r < t/ 1+ . Theorem 4 is 
able to predict the onset of such oscillations even though it says nothing about whether 
this bifurcation is supercritical or subcritical. 

Finally we perform a Turing stability analysis of the Brusselator in Section 6. We find 
that in the generic case, the modes k within the Turing instability band have the order 
0{ j) where 5 = V eD is the width of the interface. On the other hand the mesa instability 
threshold occurs when K = O ( j^r j ■ It is then clear that k K by a logarithmically 
large amount. This is the underlying reason for the coarsening process observed in Figure 
La. 

One of the easy consequences of Turing's analysis is the existence of the regime for 
which mesa-type structures are stable at the same time as the homogeneous steady state. 
A more interesting question is the following. 



Open problem 5 Does there exist a parameter regime for which the mesa-type solution 
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is unstable, and in addition the homogeneous steady state u = A, v = ^ is unstable with 
respect to the Turing instability? 

An existence of such a regime would imply spatio-temporal chaos in the Brusselator. 
We answer this question in the negative for the case when 

1 <L><C(5cxp(i j , 5 = VeD. (1.29) 



2 Steady state 

In this section we derive the asymptotics of the steady-state solution to the Brusselator 
(1.4). Let 

w = v + u. (2.1) 

and rewrite (1.4) as 

v t = 5 2 v xx + B (w - v) - (w - v) 2 v, (2.2) 
j{v t +T {w t - v t j) = Dw xx -w + v + A, (2.3) 

where 

S 2 = sD. (2.4) 
The steady-state equations then become 

f = 5 2 v + B (w-v)-(w- vf v, 

\ = Dw xx -w + v + A 1 ' J v ' 

v x {0) = = v x (L), w x (0)=0 = w x (L). (2.6) 

Next we expand (2.5) in 

v = v + ^v 1 + ..., (2.7) 

w = w + j^wi + ■ ■ ■ (2.8) 
We obtain w'q — so that w is a constant to be determined. For v we obtain 

S 2 v 0xx = F{v ,w ), (2.9) 

where 

F(v, w) = —B (v — w) + (v — w) 2 v. (2.10) 

We seek solutions for w such that F(v ,w ) is a cubic in v with equidistant roots. 
This is the so-called Maxwell line condition, and implies that the integral of F between 
its first two roots is the negative of the integral between the last two. When this is the 
case, the solution for vq will be a kink-like solution in the form of a tanh. An example 
of such solution is shown on Figure 2. The condition of equidistant roots is equivalent to 
solving two equations 

F vv (v,w) = = F(v,w) (2.11) 



Localized patterns in the Brusselator 
for unknowns B and w. A simple computation shows that this is equivalent to 



B = ^wl (2.12) 



Substituting for B we obtain 



F(v w ) = (vo - w Q f v - (v - w ) = - ^ w ^j - ^ w o^j («o - wo) ■ 

(2.13) 

On the entire space, the ODE S 2 v' ' — F(v ,wo) with F as in (2.13) admits the following 
two solutions, 

2 1 / u>o x \ . . 

v r = - W o + W o^[ Y7Ts ), (2-14) 

2 1 / wo x \ 

Vl = 3 Wo - Wo 3 t&nh {Y7Ts)- (2 ' 15) 

We are interested in mesa-type solutions. A single, symmetric mesa-type solution on an 
interval [0, L] has the form, 

^ ( v t (x - x t ) , x < % 

\ V r (x — X r ) , X > 

Here, we choose 

L ~l L + l 

xi - x r = —. (2.17) 

where I is the width of the mesa to be determined. A K-spikc symmetric solution on the 
interval [0, 1] is then obtained by glueing together K solutions on the interval L = jf. 
Such a solution has 2K interfaces whose locations are given by xu,x r i defined in (1.15). 
To find I we need the second order equations. We have, 

S 2 v lxx = F„(u ,Wo)t>i + F w (v w )w 1 , (2-18) 
w\ xx = w - v - A, (2-19) 

where 

F v (v, w) = B + (v — w) (3v — w) , F w (v, w) = —B + 2(w — v)v. (2.20) 
Integrating (2.19) and using the boundary condition w[ (±L) = we obtain 

L 

(w -v Q -A) = 0. (2.21) 



/' 

Jo 



We evaluate 

v ~l^+w (L-l) (2.22) 
from where 

, 3 LA LA ,„ . 

l~- _. 2.23 

2 wo V2B 

Substituting L = then yields Proposition 1. 

We remark that the equation (2.19) for w\ with Neumann boundary condition is 

solvable up to a constant. See Section 5, Lemma 9 for the determination of this constant. 
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3 Stability in the regime DK 2 < O (±ln 2 e) 

In this section we derive Theorem 2, valid when D <C O ( ^ ) . Before showing this 
result, we derive a more general formula for eigenvalues which is valid for all r. We show 
the following. 

Lemma 6 Consider a K-mesa symmetric equilibrium solution as given in Proposition 
1. Moreover suppose that 

1< DK 2 < O ( -ln 2 e ] . (3.1) 



The eigenvalues of such equilibrium state are asymptotically given implicitly by 



where a may take one of the following 2K values 

a j± = (c ± y/a 2 + b 2 + 2afecos (0)) , = ^ for j = 1 . . . K - 1, (3.3) 

cr ± =c + a±6, (3.4) 

w/iere 

« = ginh , b = ginh ^^ » c = /idCoth(/i d d)+/i i coth(/ijZ), (3.5) 

V2 £ + A(2t-1) VA 

W = - ;j > Md=— • (3.6) 

Proof. We start by linearizing around the equilibrium solution (v, w). We write, 

v(x,t) = v(x) + e xt (j){x) , w(x, t) = w (x) + e A V (x) (3.7) 
where we assume that ip and <f> are small. We obtain 

\<t> = 5 2 4> xx - F v (v, w)4> - F w (v, w)ip, (3.8 a) 

-{<j> + T{il>-<t>))\ = D^ xx -^ + <t>, (3.8 6) 

£ 

where 5 = \fe~D. Using (2.12) we obtain 

F v (vo,wq) = 3wq - 4w uo + (wl + B) = 3u 2 - 4u7 vo + ~jJ" w o; ( 3 - 9 ) 

2 

i^(wo, wo) = -2w 2 + 2w o«o - -B = -2w 2 + 2w Q v - -w 1 (3.10) 

so that 

F v (w ,w )=B, F w (w ,w ) = -B, (3-11) 

F v (^-,w )=B, F w (y> w °) = B - ( 3 - 12 ) 

Note that away from kink locations xu , x r i , the diffusion term eD(f>" may be neglected 

and wc have (j> ~ — f™(^°'™o) ^- < - )n ^ e other hand, near the kink locations wc locally 
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estimate the eigenfunction by the derivative of the profile. This suggests the following 
asymptotic form: 



(3.13) 



cuv u , x — xu 

ip, x (xu,x ri ) , i = l...K 
k -ip, xe(xii,x r i), i = l...K 
where the constants cu and c r i arc to be found. Wc multiply (3.8 a) by v' u and integrate. 
Because of exponential decay, we obtain that / v' H <j) ~ cu J v' t 2 . Therefore 



A 



F w (vu,w Q )tp 



1 



S v' u <t> xx - v' u <t>F v (vu,w ) - v' u 4>— (F vv (vu, w ) V! + F vw (vu,w ) wi) 



D 



(3.14) 
(3.15) 



Using integration by parts we obtain, 

J [eF)vu<j>xx ~ v' u <j>F v (vu,w )] ~ 0. 
Next we evaluate 

1 = j ' v' u <f>{F vv (v H ,w )v 1 + F vw (v H ,w )wi) . (3.16) 

Differentiating (2.18) we have 

5 2 v"-v[ (F vv (vu,w ) vi + F vw (vu,w )w 1 )-F v (v u ,w ) v[-F w (v u ,w ) w[ = (3.17) 
so that 

/ = J <f>(8 2 v'{-F v {v li ,w )v , 1 -F w (v li ,w )w' 1 ) (3.18) 
~ - J cuv'iF w (vu,w ) w[. (3.19) 



Therefore we obtain 

„+ 

Q ' 



r x u r x n r x n 

* A / v' H +ip(xu) v' H F w (vo,w ) ~ cuw'^xu) v' H F w (v ,w ) • (3.20) 

r x + 

Here and below, the symbol J L l denotes integration over the interface located at xu. 
Since v' decays exponentially outside the interface, this symbol is unambigious; that is 

f X - — fn +e.s.t. Next we show that 

x u 

J x X \' H F w (v , W0 )~-^ W l / 4 < 2 -^- (3.21) 

We have j x L v' u F w (v , w ) = G(v(xl i )) - G(v(x H )) where G(v) = J F w (v,w )dv = 
x i 

— ^v 3 +w v 2 — ^WqV. We have v{xf i ) = wo/3, v(x~^) = u>o so that G{v{xf i )) — G(v(xf i )) = 

[§T - |] W = ^^L^O- 

To evaluate the second integral in (3.21), use the explicit formula (2.15) for vi and the 
fact that J_ 00 sech 4 (y) dy = |. 
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Using (2.19) and (2.23) yields 

, w ldK 
Wi fai) g^-- (3-22) 



Therefore we obtain 



V2 / flldK 



Q»A— =4 I cu I 3-^-^0 ) + -0 (a?ii) ) , (3.23) 



An analogous computation yields 



V2 / flldK 



-c ri X— = 4 y-c ri w oJ + ip (xri)j , (3.24) 

It remains to compute V ( x u) ■ Inside the intervals (xu,x r i) we have <f) ~ an d 
outside those intervals we have <~ ip. In addition we assume that near kinks, r/> changes 
much slower than <j>. In this case we may replace 

.2 2 
% ~ -gWo<5(a; - a;;,) , u ri ~ -w S (x - x ri ) (3.25) 

inside the equation (3.8 b). Here and below, S denotes the Dirac delta function. Therefore 
we obtain, 

K 

tpxx - M 2 ^ ~ -at X! Qi< * ( x _ ^ _ Cri<5 ( x _ (3.26) 



i=l 



where 



and 



J2e + A(2r- 1) 
/x ~ w = , a; e (a; K , a; ri ) , (3.27) 

fi~/j, d = ^, x^(xu,x ri ), (3.28) 



2 (l-r)A-e ,„ „„ N 

a=3«*>- ^ • ( 3 - 29 ) 



Next we apply the following lemma. 

Lemma 7 Suppose that 

u" - [i 2 u = - (a; - x u ) + b ri S (x - x ri )j (3.30) 

with u (0) = = u (1) , (3.31) 

where xu,x r i are given in (1.15) and 

{Hi, x €E (xu,x r i) ^ 




(3.43) 
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Therefore we have 

where a are 2K eigenvalues of M -1 . 

This completes the proof of Lemma 6. ■ 

We now come back to the proof of Theorem 2. Assuming 0(t — 1) » 0, the dimensional 
analysis shows that two scalings are possible, either A = 0(e) or A = 0(De). 
We first consider the case 

A = eA . (3.45) 
Using the notation of Lemma 6 we then obtain from (3.27) and (3.28) 

rf- i±*gl -" «!. ,3.40) 



and 



i 6 ~-T' c ^ + 7 = t4' (3 ' 47) 



cr i± = c ± ^(a + 6) 2 - 2abt, t = 1 - cos (^J^j e (°> 2 ) 



(3.48) 



I // 1 \ 2 



Therefore we obtain 



m^^m)- 2 *' (3 - 49) 

^ (l ± v 7 ! - 2tf 2 d/t) . (3.50) 



Ao „ 2/^ A /A f 1 - 2 ("- 1 ) A ° + 1 ) . (3.51) 
V sD \ 1 ± Vl-2K 2 dltJ V ' 

When Aq = O (1) , the right hand side dominates and we therefore obtain 



X ° = 2(7=1) ' (3 ' 52) 

Moreover, since d = — I, I £ (0, ^) , and since i € (0, 2) we see that 2K 2 dlt < 1. This 
shows that in the case r > 1, the eigenvalues corresponding to the nodes cr,± are all real 
negative. The node cr+ is, 

o-+ = ^dtanh (^d^j + Mi tann (3.53) 

(3-54) 

~ 2^(A d+(2 + A (2r-l))Z), (3.55) 

Ao ~ 2 J*, f - ^[(, 1) A + 1] A 

Ve^V (A rf+(2 + Ao(2r-l))0/ V 7 
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Since the third term is asymptotically bigger, the assumption A = O(l) leads to 

A = — V (3.57) 

Note that this is a special case of the formula (3.52) with ± = — and t = 0. For the mode 
er_ we have 

a- = ^ d tanh (^d^j + Mi cotn (^4) ~ I' *- 3 ' 5 ^ 



Ao~2/ys^(dK--[(T-l)Ao + l]), (3.59) 

A : -■ (3.60) 

T — 1 T — 1 

Note that this is a special case of the formula (3.52) with ± = — and t = 2. 
Next we consider large eigenvalues, A 3> 0(e). Then we have 



y/(2r-l)A VA 
Mi 7= , Md~^=, (3-61) 



and we write, 

\~2y/B [ 1-y/DeldK - — - — -(t-1)Vx\. (3.62) 

\ D ) 

Dimensional analysis shows that the only way to achieve this balance is when a is large. 
But this is only possible when 

sinh (ml) = or sinh (p d d) = 0. (3.63) 

Thus either inl = inm or fi d d = inm where m is some integer. This yields the following 
eigenvalues, 

^ m 2 7r 2 , ^ m 2 7r 2 

A ~ - De WW=T) or ~ ~dP~' (3 ' 64) 

Finally, we show that a Hopf bifurcation occurs in the regime O (t — 1) <C 1. Since the 
small eigenvalues are negative when t - 1 > and positive when r — 1 <C 0, the real 
part changes sign precisely when O(t-I) « 1. To show that this occurs via a Hopf 
bifurcation, it suffices to show that A can never be zero. Suppose not. Then using some 
algebra we arrive at the following for the modes Oj± 

-1 ± y/l - 2K 2 dlt = 0, t = 1-cos (?jpj e (0,2). (3.65) 

But this is clearly impossible since K 2 dlt < \ as mentioned above. The modes a± are 
handled similarly. 

This completes the proof of Theorem 2. ■ 



4 Hopf bifurcation, DK 2 < O (±ln 2 e) 
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As Theorem 2 shows, a Hopf bifurcation occurs when r is near 1. In this section we study 
this regime in more detail, culminating in the proof of Theorem 4. 
We start by analysing the a + node. We make the assumption that 

A«0(eL>). (4.1 

This assumption will be verified later on to be consistent with the condition (1.26) 
D > \p^Q- We have, 

M ?~_<1, M3~^<1, (4-2 
and we expand cr + up to second order, 

cr + = ^ d tanh (^dT^J + Mi tann (<"4) *- 4 ' 3 

£ £>2 24 Ve£>/ eD 2 24 \ eD ' 

A \ 2 fd 3 + l 3 \ A ( 1 , 

+ -. (4.5 



eD ) \ 24 j eD \ 2K j D 
We write the equation for A as 



or 



where 



Act „ 2^ (Wa - 2 (T ' 1 £ )A + £ ) (4.6 
aA 3 + &A 2 + cA + e = 0, (4.7 

+ l ' (48 



eD J \ 24 

6 = -A f m _ 2 ay ^ , (4 . 9 

eD\2K V £> \eD \ 24 ' 



- - 



using (1.26); (4.10 



-^ B D-{^D-\^-^ 



— 4 y B ^- ( 4 - 14 

Substituting A = i\ and separating the real and imaginary part, we find that 
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This yields, 

X l = V8K(e 3 DB) 1/4 (4.16) 
and, keeping also all lower order terms for reference, 

^ * 

ID D + ^^IdK 2 (M) 

The formula (1.27) is obtained by dropping the last term of the right hand side of (4.17) 
(which is of smaller order than the first term on the right hand side), as well as dropping 
the second term in the denominator of the second term. 

Finally we must also show that the a + mode undergoes a Hopf bifurcation before all 
other modes. Let us now consider the cr_ mode. We assume here that fn ~ -^== ~ fj,^- 
and we rewrite (3.2) as 

(r A D + 1) = F(A ) = A V^o ^anh V^o ^ + coth V^o (4.18) 
where r = 1 + r , A = eD\ . (4.19) 
Near the origin and for real Ao, the curve F (Ao) ~ jAo + (f + Aq is convex and in- 
creasing. The left hand side of (4.18) is a line in A and it intersects the y axis at ~ ^/jfj 
which is a very small value by assumption (1.26). Therefore this line will intersect the 
curve F (Ao) for some small value of Ao unless its slope is precisely the slope of F(Ao) at 
the origin, i.e. — 4y^jT ~ |. This is precisely the scaling on which the Hopf bifurcation 
occurs. Subsituting I = -4= we obtain, 

Performing a similar study for the modes o~j± we obtain 



r hj± ~ 1 - J^ ARdl i 1 ± Vl - 2KHlt) , t = 1 - cos l^j G (0,2) . (4.21) 

But clearly, Th j± ,Th_ < Th + since ^-jp- 3> O (-jj) by the assumption (1.26). This shows 
that the eigenvalue A + corresponding to a + undergoes the Hopf bifurcation before any 
of the other eigenvalues, as r decreases past . ■ 

In figure 4 we show the Hopf bifurcation values and 77, _ computed numerically 
as well as the asymptotic results (1-27), (4.20), for various values of D while fixing 
S = \[T)e = 0.01. This figure shows a very good agreement when D 8. 



5 Asymmetric K-mesa solutions and instability with DK 2 <~ O (j-ln 2 e) 

In Section 3 we have shown that K mesas are always stable provided that r > 1. In our 
analysis there, we have ignored the effect of the exponentially decaying tail of v. However 
as DK 2 increases, this effect eventually must be taken into account. As we will see in this 
section, this occurs when D > O ( - ^ - ) . The main result of this section is the following. 
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~ 6 -1 -0.5 6 0^5 1 1.5 

\og 10 (D8) 

Figure 4. The value of Th + and r h _ as a function of D, while 5 = \feD = 0.01 is held fixed. The 
dots represent the numerical solution obtained by substituting A = i\i into (3.2) and solving 
for r and A; using Newton's method. The solid lines are represent formulas (1.27) and (4.20). 
Here, A = 1 and B = 15. 



Proposition 8 Let 



f(D) = - A) \3V2B (oxp + cxp {-^L= - A) } 

and Zei Di be the minimum of f(D). Let 



(5.1) 



Dk = ^I>i. (5.2) 



Suppose that r > 1 and if > 2. T/ien if — mesa solution of Proposition 1 is stable 
when D < Dk and unstable when D > Dk- The minimum D\ satisfies the following 
transcendental equation, 

(5.3) 

Suppose that 2 A 2 < B. Then 



Th8B 



A 2 

Di 7 (5-4) 

2 g ln 2 ( 12 ^ B3 [l 

e(V2B-A) 



Suppose that 2 A 2 > B. Then 

. 2 



r2B — A 

Di ^—7 ■= ' v ■ (5.5) 

2e ln 2 (i^^/ 2 ) 



Before providing a proof, consider a numerical example. Take 

s = 0.001, A = 2, B = 18. 



(5-6) 
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Then solving (5.3) we obtain Di = 21.16 so that 
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D 2 = 5.3, L> 3 = 2.35. (5.7) 

To verify Proposition 8, we ran the full numerical simulation of (1.4) for various values 
of D. We took t = 3, the initial condition to be as given in Proposition 1 with K = 2, 
and we took D from 5 to 6 with 0.1 increments every 2500 time units. For each value of 
D, we then plotted the difference of the height of the two mesas versus time. See Figure 
5. a. From this computation, we see that a change of stability occurs when D rts 5.5: if 
D < 5.5 then the difference in height is decreasing but is increasing if D > 5.5. This 
agrees well with the theoretical prediction D — 5.3. We then took K = 3, and D from 
1.9 to 3 with 0.1 increments every 2500 time units. Figure 5.b shows the the difference 
in heights of the first and second mesa. From this figure we conclude that the change in 
stability occurs when D ~ 2.45. Again, this agrees well with the theoretical prediction of 
D 3 = 2.35. 




Figure 5. (a) The difference in height of a two- mesa solution for various values of D. Here, 
D = 5 + 0.1floor(t/2500). Note the change of stability when D ~ 5.5. (b) The difference in 
height of the first two mesas of a three-mesa solution. Here, D = 1.9 + 0.1floor(t/2500). Note 
the change of stability when D ~ 2.45. In both figures, e = 0.001, A = 2, B = 18, r = 3. 



Proposition 8 follows from the existence of asymmetric patterns. Indeed a similar 
phenomena was studied for the spike solutions of the Gierer-Meinhardt model [15] and 
in the Gray-Scott model [17]. In both of these models, an asymmetric spike pattern was 
found to bifurcate from a symmetric K spike solution when K > 1. Moreover a change 
of stability of a K -spike pattern occured precisely at that point. 

To show existence of asymmetric patterns, it suffices to compute w(L) as a function 
of the domain length L, and to show an existence of a minimum of this curve. Below we 
show that such minimum occurs precisely when the the interaction in the u component is 
balanced by the interaction of v in the exponential tail. The main result is the following. 
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Lemma 9 Consider a symmetric mesa-type solution on domain [0,L]. Then we have, 

LA 1 



V2eD J 



exp 



L 



2B 



(5.8) 



w(L) 



9.04 
9.035 

9.03 
9.025 

9.02 
9.015 

9.01 
9.005 




0.4 0.5 



0.6 



0.7. 



0.8 



0.9 



Figure 6. The value of w(L) as a function of L. Here, A = 2, B = 18, e = 0.001 and D = 10. 
The solid curve represents an exact numerical value computed using the boundary value problem 
solver; the dashed curve represents the asymptotic formula (5.8). Note that both curves give 
almost the same minimum value of L w 0.68 

Proof. We recall that upon expanding the solution in as w — Wo + -j)W\ + . . ., 
v = vo + jjVi + . . ., the equation for w\ is 

S 2 v lxx = F v (v w )v 1 + F w (v w )w 1 , 



wixx =w -v - A. 



(5.9) 
(5.10) 



Note that we have 



S 2 v'oxx = F v (v ,w )v' . 



Therefore, upon multiplying (5.9) by v' and integrating we obtain, 



L/2 



v' F w (v w )wi = 5 2 {v lx v 0x - ViV 0x x)l=o /2 ■ 



To evaluate the right hand side, we write 

f-L/2 



* d 



J v' F w (v w )wi - wi (xi) J _ —(■' ; ir n ) 



(5.11) 
(5.12) 

(5.13) 



where 
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G(vo) = J F w dv = -Bv + vIwq - ^vl, (5.14) 
G(w ) = G(^) = ±wl (5.15) 



It follows that 



L/2 g 

v' F w (v w )w 1 ~ -wi (xi) gj^g. (5.16) 

To evaluate the right hand side of (5.12) we expand the solution near the boundary. 
At x = 0, we assumed that v ^ wq; writing 

v = w a + V{x) (5.17) 

we then obtain to leading order, 

5 2 V" = F v {w ,w )V + o(J^J (5.18) 
~ BV. (5.19) 
Imposing the boundary condition V (0) = we then obtain 



2 J 


f SB 1 






w +g w o exp < 




j> exp | 


{ * J 



V~K |exp|-^rj + exp| + ^rjj +0 (Jjj (5.20) 

for some constant K. To determine K, we impose the matching condition wq + V ~ 
«o + ^ui in the region -^L C i < i;. In this region we obtain 

2 Wq X-XiVB 

vq = -w - — t&n\\ — — (5.21) 

{.Pa "I 

(5.22) 

O I C I I U I 

Thus we obtain 

K = ^ w exp • ( 5 ' 23 ^ 

Moreover, for x <C x; we have 

<5 2 wi ra - B Vl - B Wl (5.24) 

-Bwi-Swi(O) (5.25) 

so that 



v 1 ~ wi (0) + DK exp | -^a; 1 . (5.26) 



22 



Theodore Kolokolnikov, Thomas Erneux and Juncheng Wei 



It follows that 



wis (0) v 0x (0)~—K 2 D, v lx (0) v 0xx (0) ~ — , (5.27) 

5 2 (wi^te ~ wiuo*x)*=o = -4B 2 Dcxp (L - o| • (5.28) 
Similarly, near x — ^ we have 

vo ~ ^ + cxp i — S cxp ^ — ^ I a; - ^ ) S (5.29) 





3 3 I <5 2 

from where we deduce 

i>i toi(0) + D^u; exp|— j^"- j exp ( ^-f- [ a; - ^ ) ) , (5.30) 

<5 2 {v\ x v 0x - WiU 0a:x ) a . =L/2 = — — exp I — I \ . (5.31) 

Therefore we obtain 

Wl (*,) ^w 3 ~ 413- D ( cxp ^ hoxp^ — Lil(l-/) ^ ) (5.32) 




(a;,) ~ 3V2B£ ( cxp \ - — j + cxp j- — (V2B - A) [> ) . (5.33) 

In the region < x < xi, we have vq ~ wo so from (2.19), we obtain 

w'l A, < x < x t . (5.34) 

It follows that 

Wl ( Xl ) - Wl (0) ~ -^xf = -± fc^) ■ (5-35) 



Using w = 3^B/2 and I = LA/V2B, we then obtain (5.8). ■ 

In Figure 6 we compare the asymptotic formula (5.8) with the numerically computed 
value for A = 2, B = 18, e = 0.001 and L> = 10. Note that the function L -> w(i) has 
a minimum at i w 0.7. This shows the existence of a fold point. Now suppose that D 
is chosen such that this minimum occurs precisely at L = -j^, with K > 1 an integer. 
Then the corresponding if -mesa equilibrium solution will have a zero eigenvalue. Since 
the exponential terms in (5.8) quickly die out as L increases, the solution becomes stable 
to the right of L = -g, and therefore unstable to the left of it. Proposition 8 is precisely 
this statement; it is obtained simply by scaling the L out and stating the existence of 
the fold point in terms of D instead. 

6 Turing analysis 

In this section we perform a Turing analysis of the homogenous steady state u = A, v = 
-j. In particular, we are interested in examining any possible connections between the 
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Turing instability regime (which leads to cosinusoidal-like patterns cos (2kirx) of mode 
k) and the localized mesa-like structures. We start by linearizing (1.4) around the steady 
state as follows, 

u = A + £e At cos (2knx) , v = — + ne xt cos {2kirx) , £,t?<1; 2k e N. (6.1) 
This yields a 2x2 eigenvalue problem for A. Its solution is given by 

A 2 -TA + A = 0, (6.2) 

where 

T = B-tA 2 -n(l + r)-e, A = r [n (n - B + A 2 ) + e {A 2 + n)] , n = 4k\ 2 eD. 

(6.3) 

Note that A„ =0 > so that the zero mode is unstable if T n=0 > or B — tA 2 > 0. 
Numerically, we observe that when the zero mode is unstable, it dominates and the 
system moves away from the equilibrium and quickly approaches a very long relaxation 
cycle, before any of the non-zero modes are activated. Therefore no spatial instability 
is observed. This leads to the following necessary condition for the Turing instability to 
appear: 

B-tA 2 < 0. (6.4) 

Provided this condition is satisfied, we have T < for all n. Therefore Turing instability 
will occur iff A < 0. In particular the second necessary condition is that 

B - A 2 > 0. (6.5) 

In this case, the most unstable mode is of the same order as the minimum of A, 

l2 B-A 2 -e ,„„, 

Shortly after the Turing instability is triggered, localized mesa-type structures appear 
due to the presence of steep gradients. In this regime, Turing instability cannot predict 
the final number K of mesas. Indeed, we have K < where 



= max \l,J-L\ (6.7) 




with D\ = O ( ln a ) as obtained in Theorem 3. It follows that as long as B — A 2 
and t - 1 > 0, we have 

K * =o (jti)> fc * = °G)' (6 - 8) 

so that K % <C fc*. Therefore we expect that shortly after the patterns appear, a coarsening 
process takes place whereby some of the resulting mesas dissapear until there are at most 
K* of them left. 

Consider an example shown on Figure l.a. Take A = 1, B = 8, e = 10~ 4 , D = 10. 
We take r = 10 to satisfy (6.4). From (5.3) we obtain D x ~ 44.5 so that = 2 and 
from (6.6) we obtain — 9. In a numerical simulation, we started from the homogenous 
steady state u = A, v = ^, perturbed by a very small random noise. A Turing instability 
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corresponding to the mode k — 7 first develops. At time t ~ 50 only six modes remain 
from which six mesas develop. One by one, these mesas are annihilated until only two 
remain, confirming our theory. We integrated the system until t = one million, but we 
do not expect any more mesa refinement since = 2. 

Formulas (6.6) together with (1.13) show that it is possible for the mesa patterns to be 
stable even as the homogenous steady state is stable. This occurs when A 2 /2 < B < A 2 
(i.e. I > \) with t > -j?. A more difficult question is whether one can find a regime in 
which both are unstable, and the system iterates between the two. Here we consider the 
case where D satisfies (1.26). In this regime the instability of a single mesa solution can 
only occur when t is near 1. Then (6.4) and (6.5) together imply that A ~ B 2 . So we 
set: 

B = A 2 + a, a<l. (6.9) 
From the condition T n=0 < we obtain 

r>l + ^ (6.10) 

and we suppose that a single mode k = 1 => n = ra* = A-k 2 sD is unstable. The 
condition A n=n „ < then leads, to leading order: 

4Tr 2 eD (eD4ir 2 - a) + eA 2 < (6.11) 

so that to leading order, 

a>^+An 2 eD (6.12) 

from where 

riT-2 ~ 1 + 0-025^ 
DAtt 2 D 

On the other hand, we have B ~ A 2 I = -k=\ and we have 

M -Ll± = 2l-2l 2 -- = V2-l-- (6.14) 
3 3 3 v ; 

so that from Theorem 4 we obtain 

r h - 1 + 0.020-^. (6.15) 
Therefore the instability of a mesa cannot follow the Turing instability since r > 77,. 



^> 1 + 7^2- 1 + 0-025-. (6.13) 



7 Discussion 

In Section 5 we were able to determine the instability thresholds without actually com- 
puting the eigenvalues; but simply by showing the existence of an asymmetric pattern 
bifurcating from a fold point. It is an open problem to find the full expression for the 
eigenvalues near this threshold. This would give a theoretical timescale for each of the 
step in the coarsening process. We expect that the unstable eigenvalue will decrease ex- 
ponentially in the distance between the mesas. This would explain the exponential time 
increase between the successive coarsening events, as observed in Figure l.a 

In Theorem 4 under condition (1.26) we have shown that as r is decreased near 1, 
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the first eigenvalue to cross the imaginary axis is A+, whose eigenfunction is even. This 
corresponds to a "breather" -type instability shown on Figure l.b. An open question is 
whether there exists a regime for which other eigenvalues undergo a Hopf bifurcation 
before the A + eigenvalue. For the Gray-Scott model it is known that a single spike can 
undergo a Hopf bifurcation due to a slow translational instability - which corresponds 
to an odd eigenfunction - whereby the center of the spike oscillates periodically [10], 
[17], [18]. The analogy of this phenomenon for a single mesa of the Brusselator would be 
the Hopf bifurcation of the A_ eigenvalue. One can also imagine spike-type solutions for 
the Brusselator simply by taking the limit I — > or cquivalcntly, B — > oo. It is an open 
problem to study this regime. 



It would be interesting to study the slow dynamics of the mesas, of which there are 
several types, corresponding to different eigenvalues. The first type is the slow transla- 
tional motion of the mesa such as seen in Figure 3 after time t ~ 2200. Similar motion 
has been analysed for the FitzHugh-Nagumo model on an infinite line [14]. In contrast 
to the Brusselator however, the FitzHugh-Nagumo model does not have a mass conser- 
vation constraint IK ~ f derived in Proposition 1, and does not undergo a coarsening 
process. In this sense the Brusselator resembles more the Cahn-Hilliard model or the 
Allen-Cahn model with mass constraint [33], [35]. However unlike the Cahn-Hilliard 
model, the Brusselator does not have a variational formulation, and the mass conserva- 
tion is only asymptotically valid. We remark that a similar phenomenon was also studied 
for Gray Scott and Gierer-Meinhardt models in the context of spike solutions [15], [9], 
[19]. 



A second type of slow instability is the mass exchange that occurs prior to mesa 
annhiliation as seen in Figure 3 at time t ~ 2000. This phenomenon also occurs in 
some flame-propagation problems [5], [34] and in the Keller-Segel model [16], where an 
exchange of mass takes place between two boundary spikes, and eventually leads to an 
annihilation of one of them. 



The coarsening process in the brusselator terminates when there are K = O ( s ^ e -i ) 
mesas left, where 5 is the characteristic width of the interface. This is in contrast to 
the the Cahn-Hilliard model, where the coarsening proceeds until all but one interface 
remains [33]. 



Localized structures far from the Turing regime are commonplace in reaction-diffusion 
systems such as the Brusselator, and provide an alternative pattern-formation mechanism 
to Turing instability. These structures appear whenever the Turing instability band is 
very large or when the diffusivity ratio of the activator and inhibitor is large. As we 
demonstrate in this work, Turing analysis cannot explain the diverse phenomena that can 
occur in this regime, such as coarsening and the "breather" -type instabilities. However 
singular perturbation tools can be successfully applied to asnwer many of these questions. 
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8 Appendix A: proof of Lemma 7 
Proof. Note that (3.30) is equivalent to solving 

u" — fifu = 0, i£U (xu, x r i) 
u"-/j, 2 d u = 0, x U (xu,x ri ) 

u'(0) = = u'(l). 
When x <G [xuXir] we have 
u = uu cosh(/x 2 (x - xu)) + Bu sinh (fi 2 {x - x u )) , x e [xux ir \ for i = 1 
where uu = u (xu) and Bu is to be found. We similarly have 
u = u ir cosh (x - x ri ))+B dl sinh (fii (x - x ri )) , a; G [x ri x l{i+1 )] for i = 

We define 

A- 1 ~ l 

a = av( i+1 ) — xu — — , 

ci = cosh , si = cosh (ml) , 
c 2 = cosh (jj, d d) , s 2 = cosh (^ d <i) . 

We have u ri — unci + Bust, from where 

u ri - u u ci . . 

Bu = tor i = 1 . . . K 

si 



and similarly = u ir Cd + B ri Sd so that 



B di = " UriCd for < = 1...JST-1. 



We also have 



b u = vl (x u ) - vl (x^) = u d (u^^Sd + B d (i-i)C d ) - hiBh 

( U U -U r u_x) c d \ U r i-UliCl 

V s d ) Si 

= — -HdU r ri-i) + (—fid + — Mi ) uu - —fnu ri for i = 2...K 

Sd \s d si J Si 

and similarly 

bri = ~— UlUu + ( — fi d + — fil ) U ri - — HdUr(i+l) for 1=1. --K -I 

si \s d si J s d 

Next note that 

u = A cosh (fix) , x € [0, xn] 
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for some constant A. Matching u (xj[) = u (x^) we then obtain 

u = U ^ X[1 ^ cosh [fix) , 



Cd/2 

i! (x n ) =u(xn)n 



s d/2 
Cd/2 
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(8.17) 
(8.18) 



where Sd/2 = sinh (jidd/2) , Cd/2 = cosh (p,dd/2) . Next we use the following identity. 
sinh(x/2) cosh (x) — 1 sinh(.x) 



to obtain 



Therefore we obtain 



cosh(x/2) sinh(a;) 1 + cosh(x) 



u (x n ) = fi Cd ^ - u (xn) . 



bn = H— -un - u' (x^) 

Sd 



Cd-1 



U r l - UnCl 



un 


i 1 




Sd 


si 


-1 


Cd 




— Hd 4 


—Vd ' 






Sd 


si 



si 



and similarly 



Kk = -—viuiK + I — m + —Hd + —m ) ui K . 



si 



si s d si 



This yields the matrix M : 



M 



a + c b 
b c a 

a c b 



where 



a = 



-fJ-d 

Sd 



b = 



si 



a c b 
b c + a 



Cd . ci 

c = —fJ-d H Hi- 

Sd si 



Consider the matrix 



Q 



a b 
b a 
a b 



b a 
a b 
b a 



(8.19) 
(8.20) 

(8.21) 
(8.22) 
(8.23) 

(8.24) 



(8.25) 



(8.26) 



(8.27) 



The eigenvalues of this matrix were computed in [29] (see Appendix B). It was found 
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that Q has the following eigenvalues, 

±^a 2 + b 2 + 2abcos(6), B =% for j = 1 . . . K - 1, (8.28) 
a + fc, a -b. (8.29) 
But we have M = (Q + c) . Therefore the eigenvalues of M are given by 

(c ± v/a 2 + & 2 + 2a6cos(6»)) , = ^ for j = 1 . . . K - 1, (8.30) 
c+a+fo, c + a-6. (8.31) 
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